#!/usr/bin/perl

# Based on: analyze.pl,v 1.9 2008/05/14 20:45:37 sherlock Exp $

# Date   : 16th October 2003
# Author: Oded Magger
# Original author : Gavin Sherlock

# License information (the MIT license)

# Copyright (c) 2003 Gavin Sherlock; Stanford University

# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files
# (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:

# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

use strict;
use warnings;
use diagnostics;

use IO::File;

use GO::TermFinder;
use GO::AnnotationProvider::AnnotationParser;
use GO::OntologyProvider::OboParser;

use GO::TermFinderReport::Text;

use GO::Utils::File    qw (GenesFromFile);
use GO::Utils::General qw (CategorizeGenes);

$|=1;

open(STDERR, ">&STDOUT");
# open(STDOUT, ">bla.txt");

###################################################################################
sub Usage{
###################################################################################

    my $message = shift;

    if (defined $message){

	print $message, "\n";

    }

    print <<USAGE;

IRRELEVANT:
This program takes a list of files, each of which contain a list of
genes, with one gene per line.  It will findTerms for the lists of
genes in each of the GO aspects, outputting the results to a file
named for the original file, but with a .terms extension.  It will only
output terms with an FDR of <= cutoff.

It will use the first supplied argument as the annotation file, the
second argument as the expected number of genes within the organism,
the third argument is the name of the obo file, and all subsequent
files as ones containing lists of genes.

Usage:

analyze.pl <annotation_file> <numGenes> <obofile> <cutoff> <clusters file>

e.g.

termFinder.pl ../t/gene_association.sgd 7200 ../t/gene_ontology_edit.obo 0.05 genes.txt

USAGE

    exit;

}

# we need at least 3 arguments, an annotation file, the number of
# genes in the genome, and a file of input genes to test

&Usage if (@ARGV < 5);

# now get our annotation file and number of genes

my $annotationFile = shift;
my $totalNum       = shift;
my $oboFile        = shift;
my $cutoffStr      = shift;
my $inputFile	   = shift;

my $cutoff		= $cutoffStr * 1.0;
my $outputFile	= $inputFile . '.out';

if ($oboFile !~ /\.obo$/){

    # require the obo file to have a .obo extension

    &Usage("Your obo file does not have a .obo extension.");

}

# now set up the objects we need

my $process   = GO::OntologyProvider::OboParser->new(ontologyFile => $oboFile,
						     aspect       => 'P');
my $component = GO::OntologyProvider::OboParser->new(ontologyFile => $oboFile,
						     aspect       => 'C');
my $function  = GO::OntologyProvider::OboParser->new(ontologyFile => $oboFile,
						     aspect       => 'F');

my $annotation = GO::AnnotationProvider::AnnotationParser->new(annotationFile=>$annotationFile);

my $termFinderP = GO::TermFinder->new(annotationProvider=> $annotation,
				      ontologyProvider  => $process,
				      totalNumGenes     => $totalNum,
				      aspect            => 'P');

my $termFinderC = GO::TermFinder->new(annotationProvider=> $annotation,
				      ontologyProvider  => $component,
				      totalNumGenes     => $totalNum,
				      aspect            => 'C');

my $termFinderF = GO::TermFinder->new(annotationProvider=> $annotation,
				      ontologyProvider  => $function,
				      totalNumGenes     => $totalNum,
				      aspect            => 'F');

# my $report = GO::TermFinderReport::Text->new();

# now go through each file

# my $report = GO::TermFinderReport::Text->new(); # TEMP TODO REMOVE

open INPUT, $inputFile or die $!;
open OUTPUT, ">$outputFile" or die $!;
my @lines = <INPUT>;
my @genes;
foreach my $line (@lines) {
	# Line is a cluster...
	@genes = split(' ', $line); # Prots in cluster.
	# print join('-', @genes); # TEMP TODO REMOVE
	# print "\n";  # TEMP TODO REMOVE

    my (@notFound, @ambiguous, @unambiguous);	

	CategorizeGenes(annotation  => $annotation,
		genes       => \@genes,
		ambiguous   => \@ambiguous,
	    unambiguous => \@unambiguous,
	    notFound    => \@notFound);

	# my $temp1 = scalar(@genes); # TEMP TODO remove
		
	@genes = (@ambiguous, @unambiguous);
	
	# my $temp2 = scalar(@genes); # TEMP TODO remove
	# print "Cluster size: $temp1, considered genes: $temp2\n"; # TEMP TODO remove
	
	# Find terms for each of the three onthologies.
	my $numHypotheses;
	my @allTerms = ();
	foreach my $termFinder ($termFinderP, $termFinderC, $termFinderF){
		my @terms = $termFinder->findTerms(genes => \@genes, calculateFDR => 1);
		@allTerms = (@allTerms, @terms);
	}
	
	# Iterate through the terms. We are looking for the term which covers the maximal amount of cluster members,
	# while the term's FDR is still below the cutoff threshold.
	# In case that several terms cover the same amount of members - takes the one which is the most statistically significant (lowest bonferroni corrected P-Value).
	# Corrected P-Val is used in this case instead of FDR (used for cutoff) because it is a more strict measurement.
	# "GO:XXXXXXX" is an 'unannotated' annotation and is thus also not considered.
	my $bestAnnot = "";
	my $maxAnnotatedMembers = 0;
	my $bestPVal = 1.0;
	
	foreach my $term (@allTerms) {
#		if (($term->{FDR_RATE} <= $cutoff) && ($term->{NODE}->goid ne "GO:XXXXXXX")) {
		if (($term->{CORRECTED_PVALUE} <= $cutoff) && ($term->{NODE}->goid ne "GO:XXXXXXX")) {
			if (($term->{NUM_ANNOTATIONS} > $maxAnnotatedMembers) ||
			    (($term->{NUM_ANNOTATIONS} == $maxAnnotatedMembers) && 
				 ($term->{CORRECTED_PVALUE} < $bestPVal))) {
				 # Found new best annotation!
				$bestAnnot = $term->{NODE}->goid;
				$maxAnnotatedMembers = $term->{NUM_ANNOTATIONS};
				$bestPVal = $term->{CORRECTED_PVALUE};
			}
		}
	}
	
	# Check if we found a ruling annotation.
	if ($maxAnnotatedMembers > 0) {
		# print "Ruling annotation: $bestAnnot\n"; # TEMP TODO remove
		print OUTPUT "$bestAnnot\n"
	}
	else {
		print OUTPUT "\n"	
	}
	
	# TEMP TODO remove
	# $numHypotheses = $report->print(pvalues  => \@allTerms,
		#   numGenes => scalar(@genes),
		 #  totalNum => $totalNum,
		 #  cutoff   => $cutoff);

}

close(INPUT);
close(OUTPUT);